########## REPLICATION: FUNCTIONS ##############

# First, I construct a function to estimate local projections based on the function by Adämmer (2019)
# in the lpirfs package, in particular the lp_lin_panel function. 
# Both functions embed the panel estimators available in the plm package.
# Adämmer's functions are much more extensive (including IVs, HP filters etc.)
# My function includes a more flexible structure for lags and leads as control variables 
# to easily estimate corrected LPs as described in the paper. 
# However, except for corrected LPs, the results can be easily replicated with the lpirfs function


########  local projections cumulative multiplier estimation ######

library(tidyverse)
library(plm)
library(lpirfs)
library(lmtest)
library(sandwich)

### extended with own state-dependency option, forward controls, (and propensity score weighting)

#' @name lp_builder
#' @title Linear Projections with State-Dependencies, Forward Controls, and 
#' @description Compute linear impulse responses with local projections by Jordà (2005).
#' @param data A data frame containing all necessary variables. Function transforms data frame into panel data set.
#' @param endog A dependent variable projected over h horizons.    
#' @param panelid Two identifiers for cross-section and time
#' @param shock A variable to shock with.
#' @param hor Number of horizons for projection.
#' @param cum_mult Estimate cumulative multipliers as y_{t+h}-y_{t-1} or regular multipliers of y{t+h}
#' @param model Panel model to estimate local projection, see vignette of plm package
#' @param effect Panel effect, see plm package
#' @param vcov_mat Select robust variance-covariance matrix workable in plm.
#' @param vcov_type Select vcov type. Can be "HC0", "HC1", "HC2"...
#' @param vcov_clus Select cluster variable for robust standard errors.
#' @param c_exog Exogenous controls with contemporaneous impact.
#' @param lag_exog Exogenous controls with lagged impact.
#' @param l_lags Number of lags.
#' @param lead_exog Exogenous controls with forward impact.
#' @param lag_aug Logical. Argument for inclusion of lagged endogenous variable. 
#' @param endog_lags Number of lags for endogenous variable.
#' @param statedep Logical. State dependent estimation of shock impact. Default is FALSE.
#' @param statevar Binary variable to indicate different states.
#' @param use_weight Logical. Specify for the use of weighting variable. Default is FALSE.
#' @param weights Supply character name of weighting variable.  There are some issues with weighting; later stage of function development

lp_builder <- function(
    data          = NULL, 
    panelid       = NULL,
    endog         = NULL,
    shock         = NULL,
    hor           = NULL,
    cum_mult      = TRUE,
    
    model         = NULL,
    effect        = NULL,
    ranmeth       = NULL,
    vcov_mat      = NULL,
    vcov_type     = NULL,
    vcov_clus     = NULL,
    
    c_exog        = NULL,
    lag_exog      = NULL,
    l_lags        = NaN, 
    lead_exog     = NULL,
    l_leads       = NaN,
    lag_aug       = FALSE,
    endog_lags    = NULL,
    
    statedep      = FALSE,
    statevar      = NULL,
    
    use_weight    = FALSE,
    wgts          = NULL
)
{
  # Check whether data_set is given
  if(is.null(data)){
    stop("You have to provide the panel data set as a pdata.frame")
  }
  
  # Check whether shock is given
  if(is.null(shock)){
    stop("You have to provide the name of the variable to shock with." )
  }
  
  # Check whether values for horizons are correct
  if(!(hor > 0) | is.nan(hor) | !(hor %% 1 == 0)){
    stop('The number of horizons has to be an integer and > 0.')
  }
  
  # Check whether panel model type is correct
  if(!model %in% c("within", "random", "ht", "between", "pooling", "fd")){
    stop("The type of the panel model has to be 'within', 'random', 'ht', 'between', 'pooling' or 'fd'. See
           the vignette of the plm package for details." )
  }
  
  # Check whether the panel effect is correctly specified
  if(!effect %in% c("individual", "time", "twoways", "nested")){
    stop("The effect introduced in the model has to be 'individual', 'time', 'twoways' or 'nested'.
           See the vignette of the plm package for details." )
  }
  
  # Check whether robust covariance estimator is correctly specified
  if(!is.null(vcov_mat)){
    if(!vcov_mat %in% c("vcov", "Vcxt", "vcovBK", "vcovDC", "vcovHC", "vcovNW", "vcovSCC")){
      stop("The choices for robust covariance estimation are 'vcov', 'vcovBK', 'vcovDC', 'vcovHC', 'vcovNW', 'vcovSCC' and 'Vcxt'.
         For details, see the vignette of the plm package and Miller (2017)." )
    }
  }
  
  # Check whether lag lengths are given if necessary
  if(!is.null(lag_exog)){
    if(is.null(l_lags)){
      stop("You have to provide the lag lengths for the exogenous data with lagged impact." )
    }
  }
  
  # Check whether lag lengths are given if necessary
  if(!is.null(lead_exog)){
    if(is.null(l_leads)){
      stop("You have to provide the lag lengths for the exogenous data with lagged impact." )
    }
  }
  
  # Check whether lag length of endogenous variable are given if necessary
  if(isTRUE(lag_aug)){
    if(is.null(endog_lags)){
      stop("You have to provide the lag lengths for the endogenous variable." )
    }
  }
  
  # Check whether lag length is numeric
  if(isTRUE(lag_aug)){
    if(!is.numeric(endog_lags)){
      stop("Lag number has to be numeric" )
    }
  }
  
  # Verify that column names are not identical
  if(length(names(data)) != length(unique(names(data)))){
    stop('Please verify that each column name is unique.')
  }
  
  # Verify that column names do not include the string pattern "lag_"
  if(length(grep("lag_", colnames(data))) >= 1){
    stop('Please do not use column names that include the string "lag_" in the name.
         This cause later naming problems')
  }
  
  
  if(isTRUE(statedep)){
    if(is.null(statevar)){
      stop("State-dependent condition must be supplied. Binary variable is necessary.")
    }
  }
  
  
  # construct formulas
  
  c_reg_names <- list()     # contemporaneous controls {t}
  c_formula <- list()
  
  lag_reg_names <- list()   # lagged controls {t-k}
  lag_formula <- list()
  
  lag_aug_reg_names <- list()  # lag augmentation  
  lag_aug_formula <- list()
  
  if(is.null(c_exog)){
    c_reg_names <- NULL
    c_formula <- NULL
  }
  
  if(is.null(lag_exog)){
    lag_reg_names <- NULL
    lag_formula <- NULL
  }
  
  if(isFALSE(lag_aug)){    
    lag_aug_formula <- NULL
    lag_aug_reg_names <- NULL   # NOTE: lag_aug is logical, because endogenous variable is already given
  }
  
  data <- plm::pdata.frame(data, index = panelid)
  
  id_vec <- data %>% dplyr::select(all_of(panelid))
  
  shock_names <- data %>% dplyr::select(all_of(shock))
  
  shock_reg_names <- colnames(shock_names) # name lag impact vector for regression formula
  
  base_formula  <- paste0(" ~ ", shock_reg_names)               
  
  if(!is.null(c_exog)){
    c_names <- data %>% dplyr::select(all_of(c_exog))
    c_reg_names <- colnames(c_names) # name contemporaneous impact vector
    c_formula   <- paste(c_reg_names, collapse = " + ")          
  } else NULL
  
  
  lag_names <- data %>% dplyr::select(all_of(panelid),all_of(lag_exog)) # extract dataframe with lagged variables
  
  
  if(!is.null(lag_exog)){
    if(ncol(lag_names)>3){
      
      ncol_vec1 <- as.numeric(ncol(lag_names))
      lag_vec   <- as.numeric(l_lags)
      name1 <- colnames(lag_names[,3:ncol_vec1])
      
      for (i in 1:lag_vec) {
        vec1 <- as.data.frame(sapply(3:ncol_vec1, function(x)plm::lag(lag_names[,x], k=i)))
        lag_names <- cbind(lag_names,vec1)
        
      }
      
      lag_names <- lag_names[,3:ncol(lag_names)]
      colnames(lag_names) <- sapply(0:l_lags, function(x)paste(name1, x, sep = "_lag"))
      lag_names <- lag_names[,(ncol_vec1-1):ncol(lag_names)]
      lag_reg_names <- colnames(lag_names) # name lag impact vector for regression formula
      lag_formula   <- paste(lag_reg_names, collapse = " + ")  
      
    } else if (ncol(lag_names)==3) {
      
      ncol_vec1 <- as.numeric(ncol(lag_names))
      lag_vec   <- as.numeric(l_lags)
      name1 <- lag_exog
      
      if(lag_vec>1){
        
        for (i in 1:lag_vec) {
          vec1 <- as.data.frame(sapply(3:ncol_vec1, function(x)plm::lag(lag_names[,x], k=i)))
          lag_names <- cbind(lag_names,vec1)
          
        }
        
        lag_names <- lag_names[,3:ncol(lag_names)]
        colnames(lag_names) <- sapply(0:l_lags, function(x)paste(name1, x, sep = "_lag"))
        lag_names <- lag_names[,(ncol_vec1-1):ncol(lag_names)]
        lag_reg_names <- colnames(lag_names) # name lag impact vector for regression formula
        lag_formula   <- paste(lag_reg_names, collapse = " + ")  
        
      } else if (lag_vec==1) {
        
        for (i in 1:lag_vec) {
          vec1 <- as.data.frame(sapply(3:ncol_vec1, function(x)plm::lag(lag_names[,x], k=i)))
          lag_names <- cbind(lag_names,vec1)
          
        }
        
        lag_names <- lag_names[,3:ncol(lag_names)]
        colnames(lag_names) <- sapply(0:l_lags, function(x)paste(name1, x, sep = "_lag"))
        lag_reg_names <- colnames(lag_names[2]) # name lag impact vector for regression formula
        help_vec <- as.vector(1:length(lag_names[,2]))
        lag_names <- cbind(lag_names,help_vec)
        lag_names <- lag_names[,2:3]
        lag_formula <- lag_reg_names 
        
      }
    }
  } 
  
  
  lag_aug_names <- data %>% dplyr::select(1,2,all_of(endog)) # extract data frame with forward variables
  
  
  if(isTRUE(lag_aug)){
    
    ncol_vec3 <- as.numeric(ncol(lag_aug_names))
    lag_aug_vec   <- as.numeric(endog_lags)
    name3 <- endog
    
    if(isTRUE(cum_mult)){
      
      vec3 <- as.data.frame(sapply(1:lag_aug_vec, function(i)plm::lag(lag_aug_names[,3], k=1)-plm::lag(lag_aug_names[,3], k=i+1)))
      lag_aug_names <- cbind(lag_aug_names,vec3)
      
    }
    
    if(isFALSE(cum_mult)){
      
      vec3 <- as.data.frame(sapply(1:lag_aug_vec, function(i)plm::lag(lag_aug_names[,3], k=i)))
      lag_aug_names <- cbind(lag_aug_names,vec3)
      
    }
    
    if(lag_aug_vec>1){
      
      lag_aug_names <- lag_aug_names[,3:ncol(lag_aug_names)]
      colnames(lag_aug_names) <- sapply(0:lag_aug_vec, function(x)paste(name3, x, sep = "_lag"))
      lag_aug_names <- lag_aug_names[,(ncol_vec3-1):ncol(lag_aug_names)]
      lag_aug_reg_names <- colnames(lag_aug_names)
      lag_aug_formula   <- paste(lag_aug_reg_names, collapse = " + ") 
      
    } 
    
    if(lag_aug_vec==1){
      
      lag_aug_names <- lag_aug_names[,3:ncol(lag_aug_names)]
      colnames(lag_aug_names) <- sapply(0:lag_aug_vec, function(x)paste(name3, x, sep = "_lag"))
      lag_aug_reg_names <- colnames(lag_aug_names[2]) # name lag impact vector for regression formula
      help_vec <- as.vector(1:length(lag_aug_names[,2]))
      lag_aug_names <- cbind(lag_aug_names,help_vec)
      lag_aug_formula <- lag_aug_reg_names 
      
    }
  } 
  
  if(isTRUE(use_weight)){
    wgtdf <- data %>% dplyr::select(all_of(wgts))
  }
  
  ### new data frame used for estimation
  
  newdat <- data.frame()
  
  if(is.null(c_exog) & is.null(lag_exog) & !isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names)} 
  if (!is.null(c_exog) & is.null(lag_exog) & !isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,c_names)}
  if (is.null(c_exog) & !is.null(lag_exog) & !isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,lag_names)}
  if (!is.null(c_exog) & !is.null(lag_exog) & !isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,c_names,lag_names)}
  
  if (is.null(c_exog) & is.null(lag_exog) & isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,lag_aug_names)}
  if (!is.null(c_exog) & is.null(lag_exog) & isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,c_names,lag_aug_names)}
  if (is.null(c_exog) & !is.null(lag_exog) & isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,lag_names,lag_aug_names)}
  if (!is.null(c_exog) & !is.null(lag_exog) & isTRUE(lag_aug) & !isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,c_names,lag_names,lag_aug_names)}
  
  if(is.null(c_exog) & is.null(lag_exog) & !isTRUE(lag_aug) & isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,wgtdf)} 
  if (!is.null(c_exog) & is.null(lag_exog) & !isTRUE(lag_aug) & isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,c_names,wgtdf)}
  if (is.null(c_exog) & !is.null(lag_exog) & !isTRUE(lag_aug) & isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,lag_names,wgtdf)}
  if (!is.null(c_exog) & !is.null(lag_exog) & !isTRUE(lag_aug) & isTRUE(use_weight)){
    newdat <- cbind(id_vec,shock_names,c_names,lag_names,wgtdf)}
  
  # lag augmentation not implemented with weights
  
  # if(!(nrow(id_vec) == nrow(shock_names) == nrow(c_names) == nrow(lag_names) == nrow(lead_names))){
  #  stop("Warning: Rows of subset dataframe differ. Check NA values of controls." )
  #}
  
  ### with state-dependency
  
  if(isTRUE(statedep)){
    
    statedf <- data %>% dplyr::select(1,2,all_of(statevar))
    
    statename1 <- colnames(statedf[3])
    st1 <- "1"
    statename0 <- colnames(statedf[3])
    st0 <- "0"
    
    statedf[,4] <- ifelse(statedf[,statevar] == 1, 1, 0)  # decompose binary state-dependent variable into two variables
    statedf[,5] <- ifelse(statedf[,statevar] == 0, 1, 0)
    
    statedf <- statedf[,4:5]
    
    colnames(statedf) <- paste(c(statename1,statename0), c(st1,st0), sep = "_")
    
    id_vec <- data %>% dplyr::select(all_of(panelid))
    
    statedf <- cbind(id_vec,statedf)
    
    statevar1_regname <- colnames(statedf[3])                # extract column names
    statevar2_regname <- colnames(statedf[4])
    
    statedf <- statedf[,3:4]
    
    newdat <- cbind(newdat,statedf)
    
    ## write state-dependent formulas with : separator for interaction syntax
    
    shock_state1 <- paste(shock, statevar1_regname, sep = ":")
    shock_state2 <- paste(shock, statevar2_regname, sep = ":")
    
    base_formula_state <- paste0(" ~ ", shock_state1, " + ", shock_state2)
    
    if(!is.null(c_exog)){
      c_formula_state1 <- paste(c_reg_names, statevar1_regname, sep = ":", collapse = " + ")
      c_formula_state2 <- paste(c_reg_names, statevar2_regname, sep = ":", collapse = " + ")
    }
    
    if(!is.null(lag_exog)){
      lag_formula_state1 <- paste(lag_reg_names, statevar1_regname, sep = ":", collapse = " + ")
      lag_formula_state2 <- paste(lag_reg_names, statevar2_regname, sep = ":", collapse = " + ")
    }
    
    if(isTRUE(lag_aug)){
      lag_aug_formula_state1 <- paste(lag_aug_reg_names, statevar1_regname, sep = ":", collapse = " + ")
      lag_aug_formula_state2 <- paste(lag_aug_reg_names, statevar2_regname, sep = ":", collapse = " + ")
    }
    
  }
  
  
  #### two conditional if statements
  
  
  # 1
  
  if(isFALSE(statedep)){
    
    if(is.null(c_reg_names) & is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula, sep = " + ")} 
    
    if(!is.null(c_reg_names) & is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula,c_formula, sep = " + ")}
    
    if(is.null(c_reg_names) & !is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula,lag_formula, sep = " + ")} 
    
    if(is.null(c_reg_names) & is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula,lag_aug_formula, sep = " + ")}
    
    if(!is.null(c_reg_names) & !is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula,c_formula,lag_formula, sep = " + ")}
    
    if(!is.null(c_reg_names) & is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula,c_formula,lag_aug_formula, sep = " + ")}
    
    if(is.null(c_reg_names) & !is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula,lag_formula,lag_aug_formula, sep = " + ")}
    
    if(!is.null(c_reg_names) & !is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula,c_formula,lag_formula,lag_aug_formula, sep = " + ")}
    
  } 
  
  
  # 2
  
  if(isTRUE(statedep)){
    if(is.null(c_reg_names) & is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state, sep = " + ")}
    
    if(!is.null(c_reg_names) & is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state,c_formula_state1,c_formula_state2, sep = " + ")}
    
    if(is.null(c_reg_names) & !is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state,lag_formula_state1,lag_formula_state2, sep = " + ")}
    
    if(is.null(c_reg_names) & is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state,lag_aug_formula_state1,lag_aug_formula_state2, sep = " + ")}
    
    if(!is.null(c_reg_names) & !is.null(lag_reg_names) & is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state,c_formula_state1,c_formula_state2,lag_formula_state1,lag_formula_state2, sep = " + ")}
    
    if(!is.null(c_reg_names) & is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state,c_formula_state1,c_formula_state2,lag_aug_formula_state1,lag_aug_formula_state2, sep = " + ")}
    
    if(is.null(c_reg_names) & !is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state,lag_formula_state1,lag_formula_state2,lag_aug_formula_state1,lag_aug_formula_state2, sep = " + ")}
    
    if(!is.null(c_reg_names) & !is.null(lag_reg_names) & !is.null(lag_aug_reg_names)){
      full_formula <- paste(base_formula_state,c_formula_state1,c_formula_state2,lag_formula_state1,lag_formula_state2,lag_aug_formula_state1,lag_aug_formula_state2, sep = " + ")}  
    
  } 
  
  
  
  # project endogenous variable over h horizons
  
  hor_fun <- hor-1 # because, horizon 1 is y_{t} - y_{t-1}
  
  endog_dat <- data %>% dplyr::select(all_of(panelid),all_of(endog))
  endog_var <- endog_dat %>% dplyr::select(all_of(endog))
  endog_name <- colnames(endog_var)
  
  endog_dat$tlag <- plm::lag(endog_dat[, 3], 1) # compute y_{t-1}
  
  # compute y_{t+h} - y_{t-1} and y_{t+h}
  
  if(isTRUE(cum_mult)){
    
    endog_proj <- as.data.frame(sapply(0:hor_fun, function(x)(plm::lead(endog_dat[,3], k=x)-endog_dat$tlag))) 
    
    colnames(endog_proj) <- sapply(0:hor_fun, function(x)paste(endog_name, x, sep = "_"))
    
  } 
  
  if(isFALSE(cum_mult)) {
    
    endog_proj <- as.data.frame(sapply(0:hor_fun, function(x)(plm::lead(endog_dat[,3], k=x))))
    
    colnames(endog_proj) <- sapply(0:hor_fun, function(x)paste(endog_name, x, sep = "_"))
    
  }
  
  
  if(nrow(newdat) == nrow(endog_proj)){
    newdat <- cbind(newdat,endog_proj)
  } else stop("Data frames not of identical length.")
  
  # extract column names of endogenous variable projected for t+h
  
  endog_reg_name <- colnames(endog_proj)
  
  # append endogenous variable to list of length h for formulas
  
  final_formula <- list()
  final_formula <- sapply(1:ncol(endog_proj), function(x)paste(endog_reg_name[x], full_formula))
  
  ## lead extension (Teulings & Zubanov 2014)
  
  lead_formula <- list()
  lead_colnames1 <- list()
  lead_colnames2 <- list()
  
  if(is.null(lead_exog)){
    lead_reg_names <- NULL
    lead_formula <- NULL
  }
  
  
  if(!is.null(lead_exog)){
    
    lead_names <- data %>% dplyr::select(1,2,all_of(lead_exog)) # extract data frame with forward variables
    
    if(ncol(lead_names)==3){
      
      lead_proj <- as.data.frame(sapply(0:hor_fun, function(x)(plm::lead(lead_names[,3], k=x))))
      namevec <- colnames(lead_names[3])
      colnames(lead_proj) <- sapply(0:hor_fun, function(x)paste(namevec, x, sep = "_"))
      
      newdat <- cbind(newdat,lead_proj)
      
      if(isFALSE(statedep)){
        
        lead_colnames <- list()
        
        for (i in 1:hor) {
          lead_colnames[i] <- colnames(lead_proj[i])
        }
        
        lead_formula <- accumulate(lead_colnames, function(a,b){paste(a,b,sep=" + ")})
        
        var1name <- lead_colnames1[1]
        
        var1name <- paste0("\\+ ", var1name)
        
        for (ii in 1:hor) {
          final_formula[ii] <- paste(final_formula[ii], lead_formula[ii], sep = " + ")
        }
        
        for (ii in 1:hor) {
          final_formula[ii] <- gsub(pattern = var1name, x = final_formula[ii], replacement = "")
        }
        
      }
      
      if(isTRUE(statedep)){
        
        lead_colnames <- list()
        lead_colnames1 <- list()
        lead_colnames2 <- list()
        
        for (i in 1:hor) {
          lead_colnames[i] <- colnames(lead_proj[i])
        }
        
        for (ii in 1:hor) {
          lead_colnames1[ii] <- paste(lead_colnames[ii], statevar1_regname, sep = ":", collapse = " + ")
          lead_colnames2[ii] <- paste(lead_colnames[ii], statevar2_regname, sep = ":", collapse = " + ")
          lead_colnames[ii] <- paste(lead_colnames1[ii], lead_colnames2[ii], sep =  " + ")
        }
        
        lead_formula <- accumulate(lead_colnames, function(a,b){paste(a,b,sep=" + ")})
        
        var1name <- lead_colnames1[1]
        var2name <- lead_colnames2[1]
        varname <- paste0("\\+ ", var1name, " \\+ ", var2name)
        
        for (ii in 1:hor) {
          final_formula[ii] <- paste(final_formula[ii], lead_formula[ii], sep = " + ")
        }
        
        for (ii in 1:hor) {
          final_formula[ii] <- gsub(pattern = varname, x = final_formula[ii], replacement = "")
        }
      }
    }
    
    if(ncol(lead_names)==4){
      
      lead_proj <- vector(mode = "list")
      colvec <- ncol(lead_names)
      
      for (i in 3:colvec) {
        lead_proj[[i]] <- as.data.frame(sapply(0:hor_fun, function(x)(plm::lead(lead_names[,i], k=x))))
      }
      
      colname1 <- rep(colnames(lead_names[3]), hor)
      colname2 <- rep(colnames(lead_names[4]), hor)
      
      colname1 <- paste(colname1, 0:hor_fun, sep = "_")
      colname2 <- paste(colname2, 0:hor_fun, sep = "_")
      
      lead_proj1 <- as.data.frame(lead_proj[[3]])
      lead_proj2 <- as.data.frame(lead_proj[[4]])
      
      colnames(lead_proj1) <- c(colname1)
      colnames(lead_proj2) <- c(colname2)
      
      newdat <- cbind(newdat,lead_proj1,lead_proj2)
      
      if(isFALSE(statedep)){
        
        lead_colnames1 <- colnames(lead_proj1)
        lead_colnames2 <- colnames(lead_proj2)
        
        var1name <- lead_colnames1[1]
        var2name <- lead_colnames2[1]
        
        lead_formula1 <- accumulate(lead_colnames1, function(a,b){paste(a,b,sep=" + ")})
        lead_formula2 <- accumulate(lead_colnames2, function(a,b){paste(a,b,sep=" + ")})
        
        var1name <- paste0("\\+ ", var1name)
        var2name <- paste0("\\+ ", var2name)
        
        for (ii in 1:hor) {
          final_formula[ii] <- paste(final_formula[ii], lead_formula1[ii], lead_formula2[ii], sep = " + ")
        }
        
        for (ii in 1:hor) {
          final_formula[ii] <- gsub(pattern = var1name, x = final_formula[ii], replacement = "")
          final_formula[ii] <- gsub(pattern = var2name, x = final_formula[ii], replacement = "")
        }
        
      }
      
      if(isTRUE(statedep)){
        
        lead_colnames1_1 <- colnames(lead_proj1)
        lead_colnames1_2 <- colnames(lead_proj1)
        lead_colnames2_1 <- colnames(lead_proj2)
        lead_colnames2_2 <- colnames(lead_proj2)
        
        lead_colnames <- list()
        
        for (ii in 1:hor) {
          lead_colnames1_1[ii] <- paste(lead_colnames1_1[ii], statevar1_regname, sep = ":", collapse = " + ")
          lead_colnames1_2[ii] <- paste(lead_colnames1_2[ii], statevar2_regname, sep = ":", collapse = " + ")
          lead_colnames2_1[ii] <- paste(lead_colnames2_1[ii], statevar1_regname, sep = ":", collapse = " + ")
          lead_colnames2_2[ii] <- paste(lead_colnames2_2[ii], statevar2_regname, sep = ":", collapse = " + ")
          
          lead_colnames[ii] <- paste(lead_colnames1_1[ii], lead_colnames1_2[ii], lead_colnames2_1[ii], lead_colnames2_2[ii], sep =  " + ")
        }
        
        lead_formula1_1 <- accumulate(lead_colnames1_1, function(a,b){paste(a,b,sep=" + ",collapse =" + ")})
        lead_formula1_2 <- accumulate(lead_colnames1_2, function(a,b){paste(a,b,sep=" + ",collapse =" + ")})
        lead_formula2_1 <- accumulate(lead_colnames2_1, function(a,b){paste(a,b,sep=" + ",collapse =" + ")})
        lead_formula2_2 <- accumulate(lead_colnames2_2, function(a,b){paste(a,b,sep=" + ",collapse =" + ")})
        
        var1_1name <- lead_colnames1_1[1]
        var1_2name <- lead_colnames1_2[1]
        var2_1name <- lead_colnames2_1[1]
        var2_2name <- lead_colnames2_2[1]
        
        for (ii in 1:hor) {
          lead_formula[ii] <- paste(lead_formula1_1[ii], lead_formula1_2[ii], lead_formula2_1[ii], lead_formula2_2[ii], sep = " + ", collapse = " + ")
        }
        
        lead_formula[1] <- paste(lead_formula1_1[1], lead_formula1_2[1], lead_formula2_1[1], lead_formula2_2[1], sep = " + ")
        
        var1_1name <- paste0("\\+ ", var1_1name)
        var1_2name <- paste0("\\+ ", var1_2name)
        var2_1name <- paste0("\\+ ", var2_1name)
        var2_2name <- paste0("\\+ ", var2_2name)
        
        for (ii in 1:hor) {
          final_formula[ii] <- paste(final_formula[ii], lead_formula[ii], sep = " + ")
        }
        
        for (ii in 1:hor) {
          final_formula[ii] <- gsub(pattern = var1_1name, x = final_formula[ii], replacement = "")
          final_formula[ii] <- gsub(pattern = var1_2name, x = final_formula[ii], replacement = "")
          final_formula[ii] <- gsub(pattern = var2_1name, x = final_formula[ii], replacement = "")
          final_formula[ii] <- gsub(pattern = var2_2name, x = final_formula[ii], replacement = "")
        }
        
      }
      
    } 
    
    if(ncol(lead_names)>4){
      stop("Lead function is restricted to two covariates so far.")
    }
    
  }
  
  
  ########  estimate with plm 
  
  newdat <- pdata.frame(newdat, index=panelid)
  
  mf <- as.vector(hor, mode = "list")
  
  if(!isTRUE(use_weight)){
    
    for (ii in (1:ncol(endog_proj))) {
      
      mf[[ii]] <- plm(formula = final_formula[[ii]],
                      data    = newdat,
                      model   = model,
                      effect  = effect,
                      random.method = ranmeth)
    } 
  }
  
  if(isTRUE(use_weight)){
    
    for (ii in (1:ncol(endog_proj))) {
      
      mf[[ii]] <- plm(formula = final_formula[[ii]],
                      data    = newdat,
                      model   = model,
                      effect  = effect,
                      random.method = ranmeth)
    }
  }
  
  
  
  ### get robust matrices
  
  robdf <- as.vector(hor, mode = "list")
  
  
  if(!is.null(vcov_type)){
    
    for (ii in 1:hor) {
      if(vcov_mat=="vcov"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcov(mf[[ii]]))}
      if (vcov_mat=="vcovBK"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovBK(mf[[ii]], type = vcov_type))}
      if (vcov_mat=="vcovDC"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovDC(mf[[ii]], type = vcov_type))}
      if (vcov_mat=="vcovHC"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovHC(mf[[ii]], type = vcov_type))}
      if (vcov_mat=="vcovNW"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovNW(mf[[ii]], type = vcov_type))}
      if (vcov_mat=="vcovSCC"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovSCC(mf[[ii]], type = vcov_type))}
      
    }
  }
  
  
  if(is.null(vcov_type)){
    
    for (ii in 1:hor) {
      if(vcov_mat=="vcov"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcov(mf[[ii]]))}
      if (vcov_mat=="vcovBK"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovBK(mf[[ii]]))}
      if (vcov_mat=="vcovDC"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovDC(mf[[ii]]))}
      if (vcov_mat=="vcovHC"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovHC(mf[[ii]]))}
      if (vcov_mat=="vcovNW"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovNW(mf[[ii]]))}
      if (vcov_mat=="vcovSCC"){
        robdf[[ii]] <- lmtest::coeftest(mf[[ii]], vcov. = vcovSCC(mf[[ii]]))}
    }
  }
  
  # summary 
  
  sumdf <- as.vector(hor, mode = "list")
  
  for (ii in 1:hor) {
    sumdf[[ii]] <- summary(mf[[ii]])
  }
  
  
  result <- list(formula     = final_formula,
                 reg_summary = sumdf,
                 reg_result  = robdf,
                 dataframe = newdat)
  
  
  class(result) <- "lp_built object"
  
  
  return(result)
}



















################################################
############# LP PLOT FUNCTION #################

# this function implements the plots used in the paper
# a model frame, the number of forecast horizons (has to match estimation) have to be provided

library(ggplot2)
library(ggpubr)

# note that some of the aesthetics (point estimates vs. ribbon) may be broken, stick to ribbon & pointest =T

plot_lp <- function(mf = NULL,
                    hor = NULL,
                    col = NULL,
                    pointest = TRUE,
                    statedep = FALSE,
                    statelabs = NULL,
                    col1 = NULL,
                    col2 = NULL,
                    lims = NULL,
                    ribbon = TRUE,
                    lty = NULL,
                    lty1 = NULL,
                    lty2 = NULL,
                    ptitle = NULL,
                    estimat = "within")
{
  
  if(isFALSE(statedep)){
    
    est <- as.vector(hor, mode = "list")
    est <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Estimate"][1])
    se <- as.vector(hor, mode = "list")
    se <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Std. Error"][1])
    horvec <- as.vector(1:hor)
    sumdf <- cbind(horvec,est,se)
    sumdf <- as.data.frame(sumdf)
    
    if(isFALSE(pointest)){
      plt <- ggplot(data = sumdf, aes(y=est,x=horvec)) + geom_line(colour=col, linewidth=1) + geom_ribbon(aes(x=horvec, ymax=est+se*1.645, ymin=est-se*1.645), fill="grey", alpha=.45) + geom_ribbon(aes(x=horvec, ymax=est+se*1.96, ymin=est-se*1.96), fill="grey", alpha=.425) + geom_hline(yintercept = 0, linetype = 2) + xlab("Horizon") + ylab("Estimate") + labs(title = "") + scale_y_continuous(limits = lims) + scale_x_continuous(breaks = c(1:hor), limits = c(1,hor)) + theme_pubclean() + theme(text = element_text(family = "Times New Roman"))
    } else if (isTRUE(pointest)) {
      plt <- ggplot(data = sumdf, aes(y=est,x=horvec)) + geom_line(colour=col, linetype=1, linewidth=0.6, alpha=1) + geom_ribbon(aes(x=horvec, ymax=est+se*1.645, ymin=est-se*1.645), fill="grey", alpha=.45) + geom_ribbon(aes(x=horvec, ymax=est+se*1.96, ymin=est-se*1.96), fill="grey", alpha=.425) + geom_line(aes(y=est), color=col, size=0.9) + geom_point(colour=col, fill=col, size=1.1, shape=21, alpha=2, stroke=1.1) + geom_hline(yintercept = 0, linetype = 2) + scale_y_continuous(limits = lims) + xlab("Horizon") + ylab("Estimate") + labs(title = ptitle) + scale_x_continuous(breaks = c(1:hor), limits = c(1,hor)) + theme_pubclean() + theme(text = element_text(family = "Times New Roman"), plot.title = element_text(hjust = 0.5, size = 11))
    }
    
  }
  
  if(isTRUE(statedep)){
    
    if(estimat=="within"){
      est1 <- as.vector(hor, mode = "list")
      est2 <- as.vector(hor, mode = "list")
      est1 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Estimate"][1])
      est2 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Estimate"][2])
      se1 <- as.vector(hor, mode = "list")
      se2 <- as.vector(hor, mode = "list")
      se1 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Std. Error"][1])
      se2 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Std. Error"][2])
    }
    
    if(estimat=="random"){
      est1 <- as.vector(hor, mode = "list")
      est2 <- as.vector(hor, mode = "list")
      est1 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Estimate"][2])
      est2 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Estimate"][3])
      se1 <- as.vector(hor, mode = "list")
      se2 <- as.vector(hor, mode = "list")
      se1 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Std. Error"][2])
      se2 <- sapply(1:hor, function(x)mf$reg_result[[x]][,"Std. Error"][3])
    }
    
    horvec <- as.vector(1:hor)
    sumdf <- cbind(horvec,est1,se1,est2,se2)
    sumdf <- as.data.frame(sumdf)
    
    colvec <- c("1"=col1,"1"=col1,"1"=col1,"1"=col1,"1"=col1,"0"=col2,"0"=col2,"0"=col2,"0"=col2,"0"=col2)  
    
    if(isFALSE(pointest)){
      colvec <- c("1"=col1,"0"=col2)  
      plt <- ggplot(data = sumdf, aes(y=est2,x=horvec)) + geom_line(color="0", linetype=1, linewidth=1.2, alpha=1) + geom_hline(yintercept = 0, linetype=2) + geom_ribbon(aes(x=horvec, ymax=est2+se2*1.645, ymin=est2-se2*1.645), fill="blue", alpha=.4) + geom_ribbon(aes(x=horvec, ymax=est2+se2*1.96, ymin=est2-se2*1.96), fill="blue", alpha=.35) + geom_line(aes(y=est2, color="0"), linewidth=1.2) + geom_line(aes(y=est1, color="1"), linewidth=1.2) + geom_ribbon(aes(x=horvec, ymax=est1+se1*1.645, ymin=est1-se1*1.645), fill="red", alpha=.4) + geom_ribbon(aes(x=horvec, ymax=est1+se1*1.96, ymin=est1-se1*1.96), fill="red", alpha=.35) + geom_line(aes(y=est1, color="1"), linewidth=1.2) + scale_y_continuous(limits = lims) + labs(color=statelabs) + xlab("Horizon") + ylab("Estimate") + labs(title = "") + scale_color_manual(values = colvec) + scale_x_continuous(breaks = c(1:hor), limits = c(1,hor)) + theme_pubclean() + theme(text = element_text(family = "Times New Roman"),legend.position = "bottom")
    } 
    if(isTRUE(pointest)) {
      if(isTRUE(ribbon)){
        plt <- ggplot(data = sumdf, aes(y=est2, x=horvec)) + geom_ribbon(aes(x=horvec, ymax=est2+se2*1.645, ymin=est2-se2*1.645), fill="blue", alpha=.3) + geom_ribbon(aes(x=horvec, ymax=est2+se2*1.96, ymin=est2-se2*1.96), fill="blue", alpha=.25) + geom_line(aes(y=est2,x=horvec), color="blue") + geom_point(aes(y=est2, x=horvec, shape="21", color="blue", fill="blue"), size=2) + geom_ribbon(aes(x=horvec, ymax=est1+se1*1.645, ymin=est1-se1*1.645), fill="red", alpha=.3) + geom_ribbon(aes(x=horvec, ymax=est1+se1*1.96, ymin=est1-se1*1.96), fill="red", alpha=.25) + geom_line(aes(y=est1,x=horvec), color="red") + geom_point(aes(y=est1, x=horvec, shape="24", color="red", fill="red"), size=2) + scale_y_continuous(limits = lims) + geom_hline(yintercept = 0, linetype = 2) + labs(shape=statelabs, color=statelabs, fill=statelabs) + scale_color_manual(labels=c("0","1"), values = c("blue","red")) + scale_shape_manual(labels=c("0","1"), values = c(21,24)) + scale_fill_manual(values = c("blue","red"), labels=c("0","1")) + xlab("Horizon") + ylab("Estimate") + labs(title = ptitle) + scale_x_continuous(breaks = c(1:hor), limits = c(1,hor)) + theme_pubclean() + theme(text = element_text(family = "Times New Roman"),legend.position = "bottom", plot.title = element_text(hjust = 0.5, size = 11))
      }
      if(isFALSE(ribbon)){
        plt <- ggplot(data = sumdf, aes(y=est1,x=horvec)) + geom_line(color="1", linetype=1, linewidth=0.6, alpha=1) + geom_line(aes(y=est1, color="1"), linewidth=0.9) + geom_point(aes(y=est1, color="1"), fill="white", size=1.1, shape=21, alpha=2, stroke=1.1) + geom_line(aes(y=est1+se1*1.96, color="1"), linetype=lty1) + geom_line(aes(y=est1-se1*1.96, color="1"), linetype=lty1) + geom_line(aes(y=est2, color="0"), linetype=1, linewidth=0.6, alpha=1) + geom_line(aes(y=est2, color="0"), linewidth=0.9) + geom_point(aes(y=est2, color="0"), fill="white", size=1.1, shape=21, alpha=2, stroke=1.1) + geom_line(aes(y=est2+se2*1.96, color="0"), linetype=lty2) + geom_line(aes(y=est2-se2*1.96, color="0"), linetype=lty2) + scale_y_continuous(limits = lims) + geom_hline(yintercept = 0, linetype = 2) + labs(color=statelabs) + xlab("Horizon") + ylab("Estimate") + labs(title = "") + scale_color_manual(values = colvec) + scale_x_continuous(breaks = c(1:hor), limits = c(1,hor)) + theme_pubclean() + theme(text = element_text(family = "Times New Roman"),legend.position = "bottom")
      }
    }
    
  }
  
  return(plt)
  
}



#########################################################################
############# IDENTIFYING RECESSIONS AND DEFICIT STATES #################


state_identifier <- function(df, var, countryname, thresh, tol){
  
  # select and filter dataset 
  
  df <- as.data.frame(df)
  
  df <- df %>% dplyr::filter(country==countryname) %>%
    dplyr::select(year,var)
  
  df <- df[complete.cases(df),]
  
  # make variable scale invariant with z-scores
  
  df$mean_var <- mean(df[,2])
  df$sd_var <- sd(df[,2])
  
  df$std_var <- (df[,var] - df$mean_var)/df$sd_var
  
  # set gamma probability space and run through dataset with logistic function
  # I generate 1000 random numbers between 0.1 and 5 to minimize gamma
  # to fulfill the condition set by the threshold
  
  gamma <- runif(1000, 0.1, 5) 
  
  hvec <- as.data.frame(sapply(1:1000, function(x)exp(-gamma[x]*df$std_var) / (1 + exp(-gamma[x]*df$std_var))))
  
  # trim observations at threshold level into binary indicators
  
  df[,6:1005] <- as.data.frame(sapply(1:1000, function(i)ifelse(hvec[,i] > thresh, 1, 0)))
  
  # clean dataset for merging
  
  df <- df %>% dplyr::select(year,starts_with("V"))
  
  df <- reshape2::melt(df, id="year")
  
  df <- df %>% group_by(variable) %>%
    dplyr::mutate(value_mean=mean(value)) %>%
    ungroup()
  
  # filter columns which fulfill threshold condition for state distinction
  
  thresh_inv = 1-thresh
  
  df <- df[df$value_mean>=(thresh_inv-tol) & df$value_mean<=(thresh_inv+tol),]
  
  df <- df %>% group_by(year) %>% dplyr::filter(value_mean==min(value_mean)) %>%
    dplyr::distinct(year,value_mean,.keep_all=T) %>%
    dplyr::select(year,variable,value) %>% rename(state=value)
  
  # prepare merging with weights on state
  
  hvec$year <- df$year
  
  hvec <- reshape2::melt(hvec, id="year")
  
  df <- left_join(df, hvec, by=c("year","variable"))
  
  df$country <- countryname
  
  df <- df %>% dplyr::select(!variable) %>%
    rename(weight=value) %>% relocate(country,.before=year)
  
  return(df)
}


######### DECOMPOSING RELATIVE ADJUSTMENT INCIDENCE ############

# TYPE A: COMPUTE SHARE OF INDIVIDUAL SPENDING COMPONENT MULTIPLIER OF TOTAL SOCIAL SPENDING MULTIPLIER
#         ESTIMATED IN BASELINE MODEL: HOW MUCH DOES EACH INDIVIDUAL SPENDING COMPONENT CONTRIBUTE TO 
#         THE TOTAL ADJUSTMENT OF SOCIAL SPENDING ACHIEVED
# TYPE B: COMPUTE PERCENTAGE SHARE OF INDIVIDUAL SPENDING COMPONENT REDUCTION OVER HORIZON H OF MEAN
#         VALUE OF SOCIAL SPENDING THROUGHOUT THE SAMPLE PERIOD. THIS DECOMPOSITION IS THEREFORE DEFINED 
#         FOR ALL POSSIBLE SPENDING COMPONENTS.


decomposer <- function(mbase, depvar, decomptype){
  
  lpb0 <- lp_builder(data = macro_lp, panelid = c("country","year"),
                     endog = depvar, shock = "SPENDING", hor = 8, cum_mult = TRUE,
                     vcov_mat = "vcovBK", model = "within", effect = "individual", 
                     c_exog = c("realgdpgr","unemp","elderly"),
                     lag_exog = c("realgdpgr","unemp"), l_lags = 2)
  
  if(decomptype == "A") {
    
    lpb0$reg_result[[1]][1] / mbase$reg_result[[1]][1] 
    lpb0$reg_result[[2]][1] / mbase$reg_result[[2]][1] 
    lpb0$reg_result[[3]][1] / mbase$reg_result[[3]][1] 
    lpb0$reg_result[[4]][1] / mbase$reg_result[[4]][1] 
    lpb0$reg_result[[5]][1] / mbase$reg_result[[5]][1] 
    lpb0$reg_result[[6]][1] / mbase$reg_result[[6]][1] 
    lpb0$reg_result[[7]][1] / mbase$reg_result[[7]][1] 
    lpb0$reg_result[[8]][1] / mbase$reg_result[[8]][1] 
    
    var_mean <- mean((macro_data[,depvar] / macro_data$socexp_public), na.rm=T)
    
    var_df <- c((lpb0$reg_result[[1]][1] / mbase$reg_result[[1]][1]),
                (lpb0$reg_result[[2]][1] / mbase$reg_result[[2]][1]),
                (lpb0$reg_result[[3]][1] / mbase$reg_result[[3]][1]),
                (lpb0$reg_result[[4]][1] / mbase$reg_result[[4]][1]),
                (lpb0$reg_result[[5]][1] / mbase$reg_result[[5]][1]),
                (lpb0$reg_result[[6]][1] / mbase$reg_result[[6]][1]),
                (lpb0$reg_result[[7]][1] / mbase$reg_result[[7]][1]),
                (lpb0$reg_result[[8]][1] / mbase$reg_result[[8]][1]))
    
    var_mean <- rep(var_mean, 8)
    
    Horizon <- c(1:8)
    
    expdata <- as.data.frame(cbind(Horizon,var_mean,var_df))
    
  } else if (decomptype == "B"){
    
    lpb0$reg_result[[1]][1] / mbase$reg_result[[1]][1] 
    lpb0$reg_result[[2]][1] / mbase$reg_result[[2]][1] 
    lpb0$reg_result[[3]][1] / mbase$reg_result[[3]][1] 
    lpb0$reg_result[[4]][1] / mbase$reg_result[[4]][1] 
    lpb0$reg_result[[5]][1] / mbase$reg_result[[5]][1] 
    lpb0$reg_result[[6]][1] / mbase$reg_result[[6]][1] 
    lpb0$reg_result[[7]][1] / mbase$reg_result[[7]][1] 
    lpb0$reg_result[[8]][1] / mbase$reg_result[[8]][1] 
    
    var_mean <- mean(macro_data[,depvar], na.rm=T)
    
    var_df <- c((lpb0$reg_result[[1]][1] / var_mean),
                (lpb0$reg_result[[2]][1] / var_mean),
                (lpb0$reg_result[[3]][1] / var_mean),
                (lpb0$reg_result[[4]][1] / var_mean),
                (lpb0$reg_result[[5]][1] / var_mean),
                (lpb0$reg_result[[6]][1] / var_mean),
                (lpb0$reg_result[[7]][1] / var_mean),
                (lpb0$reg_result[[8]][1] / var_mean))
    
    Horizon <- c(1:8)
    
    expdata <- as.data.frame(cbind(Horizon,var_df))
    
  }
  return(expdata)
}



################# further decompositions, figure iv main paper #################

## two types of decompositions: 
# a) share of baseline multiplier
# b) percentage share of average spending components 

decomposer <- function(mbase, depvar, decomptype){
  
  lpb0 <- lp_builder(data = macro_lp, panelid = c("country","year"),
                     endog = depvar, shock = "SPENDING", hor = 8, cum_mult = TRUE,
                     vcov_mat = "vcovBK", model = "within", effect = "individual", 
                     c_exog = c("realgdpgr","unemp","elderly"),
                     lag_exog = c("realgdpgr","unemp"), l_lags = 2)
  
  if(decomptype == "A") {
    
    lpb0$reg_result[[1]][1] / mbase$reg_result[[1]][1] 
    lpb0$reg_result[[2]][1] / mbase$reg_result[[2]][1] 
    lpb0$reg_result[[3]][1] / mbase$reg_result[[3]][1] 
    lpb0$reg_result[[4]][1] / mbase$reg_result[[4]][1] 
    lpb0$reg_result[[5]][1] / mbase$reg_result[[5]][1] 
    lpb0$reg_result[[6]][1] / mbase$reg_result[[6]][1] 
    lpb0$reg_result[[7]][1] / mbase$reg_result[[7]][1] 
    lpb0$reg_result[[8]][1] / mbase$reg_result[[8]][1] 
    
    var_mean <- mean((macro_data[,depvar] / macro_data$socexp_public), na.rm=T)
    
    var_df <- c((lpb0$reg_result[[1]][1] / mbase$reg_result[[1]][1]),
                (lpb0$reg_result[[2]][1] / mbase$reg_result[[2]][1]),
                (lpb0$reg_result[[3]][1] / mbase$reg_result[[3]][1]),
                (lpb0$reg_result[[4]][1] / mbase$reg_result[[4]][1]),
                (lpb0$reg_result[[5]][1] / mbase$reg_result[[5]][1]),
                (lpb0$reg_result[[6]][1] / mbase$reg_result[[6]][1]),
                (lpb0$reg_result[[7]][1] / mbase$reg_result[[7]][1]),
                (lpb0$reg_result[[8]][1] / mbase$reg_result[[8]][1]))
    
    var_mean <- rep(var_mean, 8)
    
    Horizon <- c(1:8)
    
    expdata <- as.data.frame(cbind(Horizon,var_mean,var_df))
    
  } else if (decomptype == "B"){
    
    lpb0$reg_result[[1]][1] / mbase$reg_result[[1]][1] 
    lpb0$reg_result[[2]][1] / mbase$reg_result[[2]][1] 
    lpb0$reg_result[[3]][1] / mbase$reg_result[[3]][1] 
    lpb0$reg_result[[4]][1] / mbase$reg_result[[4]][1] 
    lpb0$reg_result[[5]][1] / mbase$reg_result[[5]][1] 
    lpb0$reg_result[[6]][1] / mbase$reg_result[[6]][1] 
    lpb0$reg_result[[7]][1] / mbase$reg_result[[7]][1] 
    lpb0$reg_result[[8]][1] / mbase$reg_result[[8]][1] 
    
    var_mean <- mean(macro_data[,depvar], na.rm=T)
    
    var_df <- c((lpb0$reg_result[[1]][1] / var_mean),
                (lpb0$reg_result[[2]][1] / var_mean),
                (lpb0$reg_result[[3]][1] / var_mean),
                (lpb0$reg_result[[4]][1] / var_mean),
                (lpb0$reg_result[[5]][1] / var_mean),
                (lpb0$reg_result[[6]][1] / var_mean),
                (lpb0$reg_result[[7]][1] / var_mean),
                (lpb0$reg_result[[8]][1] / var_mean))
    
    Horizon <- c(1:8)
    
    expdata <- as.data.frame(cbind(Horizon,var_df))
    
  }
  return(expdata)
}




